%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% This program replicates Panel D of Table 3 in the paper 
%%% "Optimal Portfolio Choice with Fat Tails and Parameter
%%% Uncertainty", Journal of Financial and Quantitative Analysis
%%% (2024), Raymond Kan & Nathan Lassance 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Note: due to the use of a finite number of simulations to
%%% evaluate (k1,k2,k3), the final results for columns 't (exact)
%%% and 'Elliptical (exact)' are random and will not correspond
%%% exactly to those in the paper.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
close all

%%% Load data
load Dataset25SBTM.txt;
R = Dataset25SBTM;
[~,~,nuhat] = mlstudent(R)

%%% Define variables
nobs = 10^1;%We use nobs=10^3 to create Table 3 in the paper
gam = 1;
cbps = 10;
N = length(R(1,:));
onevec = ones(N,1);
Ttotal = length(R(:,1));

%%% Exponential integral
Expfun = @(n,x) integral(@(t)exp(-x.*t)./(t.^n),1,inf);

%%%%%%%%%%%%%%
%%% T = 60
%%%%%%%%%%%%%%

T = 60;
a1 = (T-N-2)/N;
a2 = a1*(T-N-1)*(T-N-4)/(T-2);
a3 = a2/T;
a4 = ((T-N-1)*(T-N-4)/(T*(T-2)));
rho = N/T;
M = eye(T)-(ones(T,1)*ones(1,T))./T;
NumberWindow = Ttotal-T;

%%% Initialize vectors (c for two-fund, (c1,c2) for three-fund) 
%%% Six calibrations: 1) Normal (exact), 2) t (asymp), 3) t (exact),
%%% 4) Elliptical (asymp), 5) Elliptical (exact), 6) Cross-validation
cvec1T60 = zeros(NumberWindow,1);
c1vec1T60 = zeros(NumberWindow,1);
c2vec1T60 = zeros(NumberWindow,1);
cvec2T60 = zeros(NumberWindow,1);
c1vec2T60 = zeros(NumberWindow,1);
c2vec2T60 = zeros(NumberWindow,1);
cvec3T60 = zeros(NumberWindow,1);
c1vec3T60 = zeros(NumberWindow,1);
c2vec3T60 = zeros(NumberWindow,1);
cvec4T60 = zeros(NumberWindow,1);
c1vec4T60 = zeros(NumberWindow,1);
c2vec4T60 = zeros(NumberWindow,1);
cvec5T60 = zeros(NumberWindow,1);
c1vec5T60 = zeros(NumberWindow,1);
c2vec5T60 = zeros(NumberWindow,1);
cvec6T60 = zeros(NumberWindow,1);
c1vec6T60 = zeros(NumberWindow,1);
c2vec6T60 = zeros(NumberWindow,1);
w2fvec1T60 = zeros(NumberWindow,N);
w3fvec1T60 = zeros(NumberWindow,N);
w2fvec2T60 = zeros(NumberWindow,N);
w3fvec2T60 = zeros(NumberWindow,N);
w2fvec3T60 = zeros(NumberWindow,N);
w3fvec3T60 = zeros(NumberWindow,N);
w2fvec4T60 = zeros(NumberWindow,N);
w3fvec4T60 = zeros(NumberWindow,N);
w2fvec5T60 = zeros(NumberWindow,N);
w3fvec5T60 = zeros(NumberWindow,N);
w2fvec6T60 = zeros(NumberWindow,N);
w3fvec6T60 = zeros(NumberWindow,N);
Return2fvec1T60gross = zeros(NumberWindow,1);
Return3fvec1T60gross = zeros(NumberWindow,1);
Return2fvec2T60gross = zeros(NumberWindow,1);
Return3fvec2T60gross = zeros(NumberWindow,1);
Return2fvec3T60gross = zeros(NumberWindow,1);
Return3fvec3T60gross = zeros(NumberWindow,1);
Return2fvec4T60gross = zeros(NumberWindow,1);
Return3fvec4T60gross = zeros(NumberWindow,1);
Return2fvec5T60gross = zeros(NumberWindow,1);
Return3fvec5T60gross = zeros(NumberWindow,1);
Return2fvec6T60gross = zeros(NumberWindow,1);
Return3fvec6T60gross = zeros(NumberWindow,1);
Return2fvec1T60net = zeros(NumberWindow,1);
Return3fvec1T60net = zeros(NumberWindow,1);
Return2fvec2T60net = zeros(NumberWindow,1);
Return3fvec2T60net = zeros(NumberWindow,1);
Return2fvec3T60net = zeros(NumberWindow,1);
Return3fvec3T60net = zeros(NumberWindow,1);
Return2fvec4T60net = zeros(NumberWindow,1);
Return3fvec4T60net = zeros(NumberWindow,1);
Return2fvec5T60net = zeros(NumberWindow,1);
Return3fvec5T60net = zeros(NumberWindow,1);
Return2fvec6T60net = zeros(NumberWindow,1);
Return3fvec6T60net = zeros(NumberWindow,1);

%%% Rolling window loop
for t = 1:NumberWindow

    if mod(NumberWindow-t,10) == 0
        WindowsLeftT60 = NumberWindow - t
    end
    
    %Inputs
    Xis = R(t:t+T-1,:);
    Xoos = R(t+T,:);
    muhat = mean(Xis)';
    sigmahat = cov(Xis,1);
    invsigmahat = inv(sigmahat);
    theta2hat = muhat'*invsigmahat*muhat;
    theta2hat = ((T-N-2)*theta2hat-N)/T + ((2*(theta2hat)^(N/2)...
                *((1+theta2hat)^(-(T-2)/2)))/(T*beta(N/2,(T-N)/2)*betainc(theta2hat/(1+theta2hat),N/2,(T-N)/2)));
    wghat = invsigmahat*onevec./sum(invsigmahat*onevec);
    mughat = wghat'*muhat;
    varghat = wghat'*sigmahat*wghat;
    psi2hat = muhat'*invsigmahat*muhat - mughat^2/varghat;
    psi2hat = ((T-N-1)*psi2hat-(N-1))/T + ((2*(psi2hat)^((N-1)/2)...
              *((1+psi2hat)^(-(T-2)/2)))/(T*beta((N-1)/2,(T-N+1)/2)*betainc(psi2hat/(1+psi2hat),(N-1)/2,(T-N+1)/2)));
    [~,~,nuhat] = mlstudent(Xis);
    nuhat = min(1000,nuhat);
    tauhat = sum((Xis-muhat').^2,2);
    tauhat = tauhat./mean(tauhat);

    %1) Normal (exact)
    %Two-fund
    cvec1T60(t,1) = a4*theta2hat/(theta2hat+N/T);
    w2fvec1T60(t,:) = (cvec1T60(t,1)/gam)*invsigmahat*muhat;
    Return2fvec1T60gross(t,1) = w2fvec1T60(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
            wplus(i,1) = w2fvec1T60(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w2fvec1T60(t,:)'-wplus,1);
    end 
    Return2fvec1T60net(t,1) = (1+Return2fvec1T60gross(t,1))*(1-(cbps/10000)*Turnover)-1;
    %Three-fund
    c1vec1T60(t,1) = a4*psi2hat/(psi2hat+N/T);
    c2vec1T60(t,1) = a4*(N/T)/(psi2hat+N/T);
    w3fvec1T60(t,:) = (c1vec1T60(t,1)/gam)*invsigmahat*muhat +...
                      (c2vec1T60(t,1)/gam)*mughat*invsigmahat*onevec;
    Return3fvec1T60gross(t,1) = w3fvec1T60(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
            wplus(i,1) = w3fvec1T60(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w3fvec1T60(t,:)'-wplus,1);
    end 
    Return3fvec1T60net(t,1) = (1+Return3fvec1T60gross(t,1))*(1-(cbps/10000)*Turnover)-1;

    %2) t (asymp) 
    y = @(x) (nuhat-2).*rho.*x./(2*(1-rho));
    fun = @(x) y(x).*exp(y(x)).*Expfun(nuhat/2,y(x))-rho;
    eta = fzero(fun,0.5*(1+nuhat/(nuhat-2)));
    varphi = 2*eta^2*(1-rho)/(nuhat-eta*(nuhat-2)); 
    %Two-fund
    cvec2T60(t,1) = (1-rho)^2*theta2hat/((varphi/eta)*theta2hat+rho);
    w2fvec2T60(t,:) = (cvec2T60(t,1)/gam)*invsigmahat*muhat;
    Return2fvec2T60gross(t,1) = w2fvec2T60(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
            wplus(i,1) = w2fvec2T60(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w2fvec2T60(t,:)'-wplus,1);
    end 
    Return2fvec2T60net(t,1) = (1+Return2fvec2T60gross(t,1))*(1-(cbps/10000)*Turnover)-1;
    %Three-fund
    c1vec2T60(t,1) = (1-rho)^2*psi2hat/((varphi/eta)*psi2hat+rho);
    c2vec2T60(t,1) = (1-rho)^2*(eta/varphi)*rho/((varphi/eta)*psi2hat+rho);
    w3fvec2T60(t,:) = (c1vec2T60(t,1)/gam)*invsigmahat*muhat +...
                      (c2vec2T60(t,1)/gam)*mughat*invsigmahat*onevec;
    Return3fvec2T60gross(t,1) = w3fvec2T60(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
            wplus(i,1) = w3fvec2T60(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w3fvec2T60(t,:)'-wplus,1);
    end 
    Return3fvec2T60net(t,1) = (1+Return3fvec2T60gross(t,1))*(1-(cbps/10000)*Turnover)-1;

    %3) t (exact)
    Expk1vec = zeros(nobs,1);
    Expk2vec = zeros(nobs,1);
    Expk3vec = zeros(nobs,1);
    for m = 1:nobs
        lambdahat = diag(sqrt((nuhat-2)./chi2rnd(nuhat,[T,1])));
        B = lambdahat*M*lambdahat;
        oneslambda = lambdahat*ones(T,1);
        Y = randn(T,N);
        Mat1 = inv(Y'*B*Y);
        Mat2 = Mat1^2;
        Expk1vec(m,1) = trace(Mat1);
        Expk2vec(m,1) = trace(Mat2);
        Expk3vec(m,1) = oneslambda'*Y*Mat2*Y'*oneslambda;
    end
    k1 = a1*mean(Expk1vec);
    k2 = a2*mean(Expk2vec);
    k3 = a3*mean(Expk3vec);
    %Two-fund
    cvec3T60(t,1) = a4*k1*theta2hat/(k2*theta2hat+k3*N/T);
    w2fvec3T60(t,:) = (cvec3T60(t,1)/gam)*invsigmahat*muhat;
    Return2fvec3T60gross(t,1) = w2fvec3T60(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
            wplus(i,1) = w2fvec3T60(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w2fvec3T60(t,:)'-wplus,1);
    end 
    Return2fvec3T60net(t,1) = (1+Return2fvec3T60gross(t,1))*(1-(cbps/10000)*Turnover)-1;
    %Three-fund
    c1vec3T60(t,1) = a4*k1*psi2hat/(k2*psi2hat+k3*rho);
    c2vec3T60(t,1) = a4*(k3/k2)*k1*rho/(k2*psi2hat+k3*rho);
    w3fvec3T60(t,:) = (c1vec3T60(t,1) /gam)*invsigmahat*muhat +...
                      (c2vec3T60(t,1)/gam)*mughat*invsigmahat*onevec; 
    Return3fvec3T60gross(t,1) = w3fvec3T60(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
            wplus(i,1) = w3fvec3T60(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w3fvec3T60(t,:)'-wplus,1);
    end 
    Return3fvec3T60net(t,1) = (1+Return3fvec3T60gross(t,1))*(1-(cbps/10000)*Turnover)-1;

    %4) Elliptical (asymp)
    fun = @(x) sum(1./(T-N+x.*tauhat.*N))-1;
    eta = fzero(fun,[1,10]);
    varphi = (1-rho)./(1/eta^2-sum((N.*tauhat.^2)./(T-N+N.*eta.*tauhat).^2));
    %Two-fund
    cvec4T60(t,1) = (1-rho)^2*theta2hat/((varphi/eta)*theta2hat+rho);
    w2fvec4T60(t,:) = (cvec4T60(t,1)/gam)*invsigmahat*muhat;
    Return2fvec4T60gross(t,1) = w2fvec4T60(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
            wplus(i,1) = w2fvec4T60(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w2fvec4T60(t,:)'-wplus,1);
    end 
    Return2fvec4T60net(t,1) = (1+Return2fvec4T60gross(t,1))*(1-(cbps/10000)*Turnover)-1;
    %Three-fund
    c1vec4T60(t,1) = (1-rho)^2*psi2hat/((varphi/eta)*psi2hat+rho);
    c2vec4T60(t,1) = (1-rho)^2*(eta/varphi)*rho/((varphi/eta)*psi2hat+rho);
    w3fvec4T60(t,:) = (c1vec4T60(t,1)/gam)*invsigmahat*muhat +...
                      (c2vec4T60(t,1)/gam)*mughat*invsigmahat*onevec;
    Return3fvec4T60gross(t,1) = w3fvec4T60(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
            wplus(i,1) = w3fvec4T60(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w3fvec4T60(t,:)'-wplus,1);
    end 
    Return3fvec4T60net(t,1) = (1+Return3fvec4T60gross(t,1))*(1-(cbps/10000)*Turnover)-1;
    
    %5) Elliptical (exact)
    lambdahat = diag(sqrt(tauhat));
    B = lambdahat*M*lambdahat;
    oneslambda = lambdahat*ones(T,1);
    Expk1vec = zeros(nobs,1);
    Expk2vec = zeros(nobs,1);
    Expk3vec = zeros(nobs,1);
    for m = 1:nobs
        Y = randn(T,N);
        Mat1 = inv(Y'*B*Y);
        Mat2 = Mat1^2;
        Expk1vec(m,1) = trace(Mat1);
        Expk2vec(m,1) = trace(Mat2);
        Expk3vec(m,1) = oneslambda'*Y*Mat2*Y'*oneslambda;
    end
    k1 = a1*mean(Expk1vec);
    k2 = a2*mean(Expk2vec);
    k3 = a3*mean(Expk3vec);
    %Two-fund
    cvec5T60(t,1) = a4*k1*theta2hat/(k2*theta2hat+k3*N/T);
    w2fvec5T60(t,:) = (cvec5T60(t,1)/gam)*invsigmahat*muhat;
    Return2fvec5T60gross(t,1) = w2fvec5T60(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
            wplus(i,1) = w2fvec5T60(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w2fvec5T60(t,:)'-wplus,1);
    end 
    Return2fvec5T60net(t,1) = (1+Return2fvec5T60gross(t,1))*(1-(cbps/10000)*Turnover)-1;
    %Three-fund
    c1vec5T60(t,1)  = a4*k1*psi2hat/(k2*psi2hat+k3*rho);
    c2vec5T60(t,1) = a4*(k3/k2)*k1*rho/(k2*psi2hat+k3*rho);
    w3fvec5T60(t,:) = (c1vec5T60(t,1) /gam)*invsigmahat*muhat +...
                      (c2vec5T60(t,1)/gam)*mughat*invsigmahat*onevec; 
    Return3fvec5T60gross(t,1) = w3fvec5T60(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
            wplus(i,1) = w3fvec5T60(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w3fvec5T60(t,:)'-wplus,1);
    end 
    Return3fvec5T60net(t,1) = (1+Return3fvec5T60gross(t,1))*(1-(cbps/10000)*Turnover)-1;

    %6) Cross-validation 
    kfold = 5;
    %Two-fund
    dc = 0.01;
    ckvec = 0:dc:1;
    Uk = zeros(length(ckvec),1);
    i = 1;
    for c = ckvec   
        Pk = zeros(T,1);
        for k = 1:kfold
            Xk = Xis;
            Xk(1+(k-1)*T/kfold:k*T/kfold,:) = [];
            Xoosk = Xis(1+(k-1)*T/kfold:k*T/kfold,:);
            muk = mean(Xk)';
            sigmak = cov(Xk,1);
            wk = (c/gam)*inv(sigmak)*muk;
            Pk(1+(k-1)*T/kfold:k*T/kfold,1) = Xoosk*wk;
         end
    Uk(i,1) = mean(Pk) - (gam/2)*var(Pk);
    i = i+1;
    end
    [~,index] = max(Uk);
    cvec6T60(t,1) = (index-1)*dc;
    w2fvec6T60(t,:) = (cvec6T60(t,1)/gam)*invsigmahat*muhat;
    Return2fvec6T60gross(t,1) = w2fvec6T60(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
            wplus(i,1) = w2fvec6T60(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w2fvec6T60(t,:)'-wplus,1);
    end 
    Return2fvec6T60net(t,1) = (1+Return2fvec6T60gross(t,1))*(1-(cbps/10000)*Turnover)-1;
    %Three-fund
    dc = 0.025;
    c1kvec = 0:dc:1;
    c2kvec = 0:dc:1;
    Uk = zeros(length(c1kvec),length(c2kvec));
    i = 1;
    j = 1;
    for c1 = c1kvec   
        for c2 = c2kvec
            Pk = zeros(T,1);
            for k = 1:kfold
                Xk = Xis;
                Xk(1+(k-1)*T/kfold:k*T/kfold,:) = [];
                Xoosk = Xis(1+(k-1)*T/kfold:k*T/kfold,:);
                muk = mean(Xk)';
                sigmak = cov(Xk,1);
                wgk = inv(sigmak)*onevec./sum(inv(sigmak)*onevec);
                mugk = wgk'*muk;
                wk = (c1/gam)*inv(sigmak)*muk + (c2/gam)*mugk*inv(sigmak)*onevec;
                Pk(1+(k-1)*T/kfold:k*T/kfold,1) = Xoosk*wk;
            end
            Uk(i,j) = mean(Pk) - (gam/2)*var(Pk);
            j = j+1;
        end
        j = 1;
        i = i+1;
    end
    [maxValue,~] = max(Uk(:));
    [index1,index2] = find(Uk == maxValue);
    c1vec6T60(t,1) = c1kvec(1) + (index1-1)*dc;
    c2vec6T60(t,1) = c2kvec(1) + (index2-1)*dc;
    w3fvec6T60(t,:) = (c1vec6T60(t,1)/gam)*invsigmahat*muhat + ...
                      (c2vec6T60(t,1)/gam)*mughat*invsigmahat*onevec;
    Return3fvec6T60gross(t,1) = w3fvec6T60(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
             wplus(i,1) = w3fvec6T60(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w3fvec6T60(t,:)'-wplus,1);
    end 
    Return3fvec6T60net(t,1) = (1+Return3fvec6T60gross(t,1))*(1-(cbps/10000)*Turnover)-1;

end

EU2f1T60gross = 12*(mean(Return2fvec1T60gross)-(gam/2)*var(Return2fvec1T60gross));
EU3f1T60gross = 12*(mean(Return3fvec1T60gross)-(gam/2)*var(Return3fvec1T60gross));
EU2f2T60gross = 12*(mean(Return2fvec2T60gross)-(gam/2)*var(Return2fvec2T60gross));
EU3f2T60gross = 12*(mean(Return3fvec2T60gross)-(gam/2)*var(Return3fvec2T60gross));
EU2f3T60gross = 12*(mean(Return2fvec3T60gross)-(gam/2)*var(Return2fvec3T60gross));
EU3f3T60gross = 12*(mean(Return3fvec3T60gross)-(gam/2)*var(Return3fvec3T60gross));
EU2f4T60gross = 12*(mean(Return2fvec4T60gross)-(gam/2)*var(Return2fvec4T60gross));
EU3f4T60gross = 12*(mean(Return3fvec4T60gross)-(gam/2)*var(Return3fvec4T60gross));
EU2f5T60gross = 12*(mean(Return2fvec5T60gross)-(gam/2)*var(Return2fvec5T60gross));
EU3f5T60gross = 12*(mean(Return3fvec5T60gross)-(gam/2)*var(Return3fvec5T60gross));
EU2f6T60gross = 12*(mean(Return2fvec6T60gross)-(gam/2)*var(Return2fvec6T60gross));
EU3f6T60gross = 12*(mean(Return3fvec6T60gross)-(gam/2)*var(Return3fvec6T60gross));
EU2f1T60net = 12*(mean(Return2fvec1T60net)-(gam/2)*var(Return2fvec1T60net));
EU3f1T60net = 12*(mean(Return3fvec1T60net)-(gam/2)*var(Return3fvec1T60net));
EU2f2T60net = 12*(mean(Return2fvec2T60net)-(gam/2)*var(Return2fvec2T60net));
EU3f2T60net = 12*(mean(Return3fvec2T60net)-(gam/2)*var(Return3fvec2T60net));
EU2f3T60net = 12*(mean(Return2fvec3T60net)-(gam/2)*var(Return2fvec3T60net));
EU3f3T60net = 12*(mean(Return3fvec3T60net)-(gam/2)*var(Return3fvec3T60net));
EU2f4T60net = 12*(mean(Return2fvec4T60net)-(gam/2)*var(Return2fvec4T60net));
EU3f4T60net = 12*(mean(Return3fvec4T60net)-(gam/2)*var(Return3fvec4T60net));
EU2f5T60net = 12*(mean(Return2fvec5T60net)-(gam/2)*var(Return2fvec5T60net));
EU3f5T60net = 12*(mean(Return3fvec5T60net)-(gam/2)*var(Return3fvec5T60net));
EU2f6T60net = 12*(mean(Return2fvec6T60net)-(gam/2)*var(Return2fvec6T60net));
EU3f6T60net = 12*(mean(Return3fvec6T60net)-(gam/2)*var(Return3fvec6T60net));

Table3PanelDTwoFundT60 = ...
[EU2f1T60gross,EU2f2T60gross,EU2f3T60gross,EU2f4T60gross,EU2f5T60gross,...
EU2f6T60gross;EU2f1T60net,EU2f2T60net,EU2f3T60net,EU2f4T60net,EU2f5T60net,...
EU2f6T60net;mean(cvec1T60),mean(cvec2T60),mean(cvec3T60),mean(cvec4T60),...
mean(cvec5T60),mean(cvec6T60)]
Table3PanelDThreeFundT60 = ...
[EU3f1T60gross,EU3f2T60gross,EU3f3T60gross,EU3f4T60gross,EU3f5T60gross,...
EU3f6T60gross;EU3f1T60net,EU3f2T60net,EU3f3T60net,EU3f4T60net,EU3f5T60net,...
EU3f6T60net;mean(c1vec1T60),mean(c1vec2T60),mean(c1vec3T60),mean(c1vec4T60),...
mean(c1vec5T60),mean(c1vec6T60);mean(c2vec1T60),mean(c2vec2T60),mean(c2vec3T60),...
mean(c2vec4T60),mean(c2vec5T60),mean(c2vec6T60)]

%%%%%%%%%%%%%%
%%% T = 120
%%%%%%%%%%%%%%

T = 120;
a1 = (T-N-2)/N;
a2 = a1*(T-N-1)*(T-N-4)/(T-2);
a3 = a2/T;
a4 = ((T-N-1)*(T-N-4)/(T*(T-2)));
rho = N/T;
M = eye(T)-(ones(T,1)*ones(1,T))./T;
NumberWindow = Ttotal-T;

%%% Initialize vectors (c for two-fund, (c1,c2) for three-fund) 
%%% Six calibrations: 1) Normal (exact), 2) t (asymp), 3) t (exact),
%%% 4) Elliptical (asymp), 5) Elliptical (exact), 6) Cross-validation
cvec1T120 = zeros(NumberWindow,1);
c1vec1T120 = zeros(NumberWindow,1);
c2vec1T120 = zeros(NumberWindow,1);
cvec2T120 = zeros(NumberWindow,1);
c1vec2T120 = zeros(NumberWindow,1);
c2vec2T120 = zeros(NumberWindow,1);
cvec3T120 = zeros(NumberWindow,1);
c1vec3T120 = zeros(NumberWindow,1);
c2vec3T120 = zeros(NumberWindow,1);
cvec4T120 = zeros(NumberWindow,1);
c1vec4T120 = zeros(NumberWindow,1);
c2vec4T120 = zeros(NumberWindow,1);
cvec5T120 = zeros(NumberWindow,1);
c1vec5T120 = zeros(NumberWindow,1);
c2vec5T120 = zeros(NumberWindow,1);
cvec6T120 = zeros(NumberWindow,1);
c1vec6T120 = zeros(NumberWindow,1);
c2vec6T120 = zeros(NumberWindow,1);
w2fvec1T120 = zeros(NumberWindow,N);
w3fvec1T120 = zeros(NumberWindow,N);
w2fvec2T120 = zeros(NumberWindow,N);
w3fvec2T120 = zeros(NumberWindow,N);
w2fvec3T120 = zeros(NumberWindow,N);
w3fvec3T120 = zeros(NumberWindow,N);
w2fvec4T120 = zeros(NumberWindow,N);
w3fvec4T120 = zeros(NumberWindow,N);
w2fvec5T120 = zeros(NumberWindow,N);
w3fvec5T120 = zeros(NumberWindow,N);
w2fvec6T120 = zeros(NumberWindow,N);
w3fvec6T120 = zeros(NumberWindow,N);
Return2fvec1T120gross = zeros(NumberWindow,1);
Return3fvec1T120gross = zeros(NumberWindow,1);
Return2fvec2T120gross = zeros(NumberWindow,1);
Return3fvec2T120gross = zeros(NumberWindow,1);
Return2fvec3T120gross = zeros(NumberWindow,1);
Return3fvec3T120gross = zeros(NumberWindow,1);
Return2fvec4T120gross = zeros(NumberWindow,1);
Return3fvec4T120gross = zeros(NumberWindow,1);
Return2fvec5T120gross = zeros(NumberWindow,1);
Return3fvec5T120gross = zeros(NumberWindow,1);
Return2fvec6T120gross = zeros(NumberWindow,1);
Return3fvec6T120gross = zeros(NumberWindow,1);
Return2fvec1T120net = zeros(NumberWindow,1);
Return3fvec1T120net = zeros(NumberWindow,1);
Return2fvec2T120net = zeros(NumberWindow,1);
Return3fvec2T120net = zeros(NumberWindow,1);
Return2fvec3T120net = zeros(NumberWindow,1);
Return3fvec3T120net = zeros(NumberWindow,1);
Return2fvec4T120net = zeros(NumberWindow,1);
Return3fvec4T120net = zeros(NumberWindow,1);
Return2fvec5T120net = zeros(NumberWindow,1);
Return3fvec5T120net = zeros(NumberWindow,1);
Return2fvec6T120net = zeros(NumberWindow,1);
Return3fvec6T120net = zeros(NumberWindow,1);

%%% Rolling window loop
for t = 1:NumberWindow

    if mod(NumberWindow-t,10) == 0
        WindowsLeftT120 = NumberWindow - t
    end
    
    %Inputs
    Xis = R(t:t+T-1,:);
    Xoos = R(t+T,:);
    muhat = mean(Xis)';
    sigmahat = cov(Xis,1);
    invsigmahat = inv(sigmahat);
    theta2hat = muhat'*invsigmahat*muhat;
    theta2hat = ((T-N-2)*theta2hat-N)/T + ((2*(theta2hat)^(N/2)...
                *((1+theta2hat)^(-(T-2)/2)))/(T*beta(N/2,(T-N)/2)*betainc(theta2hat/(1+theta2hat),N/2,(T-N)/2)));
    wghat = invsigmahat*onevec./sum(invsigmahat*onevec);
    mughat = wghat'*muhat;
    varghat = wghat'*sigmahat*wghat;
    psi2hat = muhat'*invsigmahat*muhat - mughat^2/varghat;
    psi2hat = ((T-N-1)*psi2hat-(N-1))/T + ((2*(psi2hat)^((N-1)/2)...
              *((1+psi2hat)^(-(T-2)/2)))/(T*beta((N-1)/2,(T-N+1)/2)*betainc(psi2hat/(1+psi2hat),(N-1)/2,(T-N+1)/2)));
    [~,~,nuhat] = mlstudent(Xis);
    nuhat = min(1000,nuhat);
    tauhat = sum((Xis-muhat').^2,2);
    tauhat = tauhat./mean(tauhat);

    %1) Normal (exact)
    %Two-fund
    cvec1T120(t,1) = a4*theta2hat/(theta2hat+N/T);
    w2fvec1T120(t,:) = (cvec1T120(t,1)/gam)*invsigmahat*muhat;
    Return2fvec1T120gross(t,1) = w2fvec1T120(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
            wplus(i,1) = w2fvec1T120(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w2fvec1T120(t,:)'-wplus,1);
    end 
    Return2fvec1T120net(t,1) = (1+Return2fvec1T120gross(t,1))*(1-(cbps/10000)*Turnover)-1;
    %Three-fund
    c1vec1T120(t,1) = a4*psi2hat/(psi2hat+N/T);
    c2vec1T120(t,1) = a4*(N/T)/(psi2hat+N/T);
    w3fvec1T120(t,:) = (c1vec1T120(t,1)/gam)*invsigmahat*muhat +...
                      (c2vec1T120(t,1)/gam)*mughat*invsigmahat*onevec;
    Return3fvec1T120gross(t,1) = w3fvec1T120(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
            wplus(i,1) = w3fvec1T120(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w3fvec1T120(t,:)'-wplus,1);
    end 
    Return3fvec1T120net(t,1) = (1+Return3fvec1T120gross(t,1))*(1-(cbps/10000)*Turnover)-1;

    %2) t (asymp) 
    y = @(x) (nuhat-2).*rho.*x./(2*(1-rho));
    fun = @(x) y(x).*exp(y(x)).*Expfun(nuhat/2,y(x))-rho;
    eta = fzero(fun,0.5*(1+nuhat/(nuhat-2)));
    varphi = 2*eta^2*(1-rho)/(nuhat-eta*(nuhat-2)); 
    %Two-fund
    cvec2T120(t,1) = (1-rho)^2*theta2hat/((varphi/eta)*theta2hat+rho);
    w2fvec2T120(t,:) = (cvec2T120(t,1)/gam)*invsigmahat*muhat;
    Return2fvec2T120gross(t,1) = w2fvec2T120(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
            wplus(i,1) = w2fvec2T120(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w2fvec2T120(t,:)'-wplus,1);
    end 
    Return2fvec2T120net(t,1) = (1+Return2fvec2T120gross(t,1))*(1-(cbps/10000)*Turnover)-1;
    %Three-fund
    c1vec2T120(t,1) = (1-rho)^2*psi2hat/((varphi/eta)*psi2hat+rho);
    c2vec2T120(t,1) = (1-rho)^2*(eta/varphi)*rho/((varphi/eta)*psi2hat+rho);
    w3fvec2T120(t,:) = (c1vec2T120(t,1)/gam)*invsigmahat*muhat +...
                      (c2vec2T120(t,1)/gam)*mughat*invsigmahat*onevec;
    Return3fvec2T120gross(t,1) = w3fvec2T120(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
            wplus(i,1) = w3fvec2T120(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w3fvec2T120(t,:)'-wplus,1);
    end 
    Return3fvec2T120net(t,1) = (1+Return3fvec2T120gross(t,1))*(1-(cbps/10000)*Turnover)-1;

    %3) t (exact)
    Expk1vec = zeros(nobs,1);
    Expk2vec = zeros(nobs,1);
    Expk3vec = zeros(nobs,1);
    for m = 1:nobs
        lambdahat = diag(sqrt((nuhat-2)./chi2rnd(nuhat,[T,1])));
        B = lambdahat*M*lambdahat;
        oneslambda = lambdahat*ones(T,1);
        Y = randn(T,N);
        Mat1 = inv(Y'*B*Y);
        Mat2 = Mat1^2;
        Expk1vec(m,1) = trace(Mat1);
        Expk2vec(m,1) = trace(Mat2);
        Expk3vec(m,1) = oneslambda'*Y*Mat2*Y'*oneslambda;
    end
    k1 = a1*mean(Expk1vec);
    k2 = a2*mean(Expk2vec);
    k3 = a3*mean(Expk3vec);
    %Two-fund
    cvec3T120(t,1) = a4*k1*theta2hat/(k2*theta2hat+k3*N/T);
    w2fvec3T120(t,:) = (cvec3T120(t,1)/gam)*invsigmahat*muhat;
    Return2fvec3T120gross(t,1) = w2fvec3T120(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
            wplus(i,1) = w2fvec3T120(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w2fvec3T120(t,:)'-wplus,1);
    end 
    Return2fvec3T120net(t,1) = (1+Return2fvec3T120gross(t,1))*(1-(cbps/10000)*Turnover)-1;
    %Three-fund
    c1vec3T120(t,1) = a4*k1*psi2hat/(k2*psi2hat+k3*rho);
    c2vec3T120(t,1) = a4*(k3/k2)*k1*rho/(k2*psi2hat+k3*rho);
    w3fvec3T120(t,:) = (c1vec3T120(t,1) /gam)*invsigmahat*muhat +...
                      (c2vec3T120(t,1)/gam)*mughat*invsigmahat*onevec; 
    Return3fvec3T120gross(t,1) = w3fvec3T120(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
            wplus(i,1) = w3fvec3T120(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w3fvec3T120(t,:)'-wplus,1);
    end 
    Return3fvec3T120net(t,1) = (1+Return3fvec3T120gross(t,1))*(1-(cbps/10000)*Turnover)-1;

    %4) Elliptical (asymp)
    fun = @(x) sum(1./(T-N+x.*tauhat.*N))-1;
    eta = fzero(fun,[1,10]);
    varphi = (1-rho)./(1/eta^2-sum((N.*tauhat.^2)./(T-N+N.*eta.*tauhat).^2));
    %Two-fund
    cvec4T120(t,1) = (1-rho)^2*theta2hat/((varphi/eta)*theta2hat+rho);
    w2fvec4T120(t,:) = (cvec4T120(t,1)/gam)*invsigmahat*muhat;
    Return2fvec4T120gross(t,1) = w2fvec4T120(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
            wplus(i,1) = w2fvec4T120(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w2fvec4T120(t,:)'-wplus,1);
    end 
    Return2fvec4T120net(t,1) = (1+Return2fvec4T120gross(t,1))*(1-(cbps/10000)*Turnover)-1;
    %Three-fund
    c1vec4T120(t,1) = (1-rho)^2*psi2hat/((varphi/eta)*psi2hat+rho);
    c2vec4T120(t,1) = (1-rho)^2*(eta/varphi)*rho/((varphi/eta)*psi2hat+rho);
    w3fvec4T120(t,:) = (c1vec4T120(t,1)/gam)*invsigmahat*muhat +...
                      (c2vec4T120(t,1)/gam)*mughat*invsigmahat*onevec;
    Return3fvec4T120gross(t,1) = w3fvec4T120(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
            wplus(i,1) = w3fvec4T120(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w3fvec4T120(t,:)'-wplus,1);
    end 
    Return3fvec4T120net(t,1) = (1+Return3fvec4T120gross(t,1))*(1-(cbps/10000)*Turnover)-1;
    
    %5) Elliptical (exact)
    lambdahat = diag(sqrt(tauhat));
    B = lambdahat*M*lambdahat;
    oneslambda = lambdahat*ones(T,1);
    Expk1vec = zeros(nobs,1);
    Expk2vec = zeros(nobs,1);
    Expk3vec = zeros(nobs,1);
    for m = 1:nobs
        Y = randn(T,N);
        Mat1 = inv(Y'*B*Y);
        Mat2 = Mat1^2;
        Expk1vec(m,1) = trace(Mat1);
        Expk2vec(m,1) = trace(Mat2);
        Expk3vec(m,1) = oneslambda'*Y*Mat2*Y'*oneslambda;
    end
    k1 = a1*mean(Expk1vec);
    k2 = a2*mean(Expk2vec);
    k3 = a3*mean(Expk3vec);
    %Two-fund
    cvec5T120(t,1) = a4*k1*theta2hat/(k2*theta2hat+k3*N/T);
    w2fvec5T120(t,:) = (cvec5T120(t,1)/gam)*invsigmahat*muhat;
    Return2fvec5T120gross(t,1) = w2fvec5T120(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
            wplus(i,1) = w2fvec5T120(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w2fvec5T120(t,:)'-wplus,1);
    end 
    Return2fvec5T120net(t,1) = (1+Return2fvec5T120gross(t,1))*(1-(cbps/10000)*Turnover)-1;
    %Three-fund
    c1vec5T120(t,1)  = a4*k1*psi2hat/(k2*psi2hat+k3*rho);
    c2vec5T120(t,1) = a4*(k3/k2)*k1*rho/(k2*psi2hat+k3*rho);
    w3fvec5T120(t,:) = (c1vec5T120(t,1) /gam)*invsigmahat*muhat +...
                      (c2vec5T120(t,1)/gam)*mughat*invsigmahat*onevec; 
    Return3fvec5T120gross(t,1) = w3fvec5T120(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
            wplus(i,1) = w3fvec5T120(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w3fvec5T120(t,:)'-wplus,1);
    end 
    Return3fvec5T120net(t,1) = (1+Return3fvec5T120gross(t,1))*(1-(cbps/10000)*Turnover)-1;

    %6) Cross-validation 
    kfold = 5;
    %Two-fund
    dc = 0.01;
    ckvec = 0:dc:1;
    Uk = zeros(length(ckvec),1);
    i = 1;
    for c = ckvec   
        Pk = zeros(T,1);
        for k = 1:kfold
            Xk = Xis;
            Xk(1+(k-1)*T/kfold:k*T/kfold,:) = [];
            Xoosk = Xis(1+(k-1)*T/kfold:k*T/kfold,:);
            muk = mean(Xk)';
            sigmak = cov(Xk,1);
            wk = (c/gam)*inv(sigmak)*muk;
            Pk(1+(k-1)*T/kfold:k*T/kfold,1) = Xoosk*wk;
         end
    Uk(i,1) = mean(Pk) - (gam/2)*var(Pk);
    i = i+1;
    end
    [~,index] = max(Uk);
    cvec6T120(t,1) = (index-1)*dc;
    w2fvec6T120(t,:) = (cvec6T120(t,1)/gam)*invsigmahat*muhat;
    Return2fvec6T120gross(t,1) = w2fvec6T120(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
            wplus(i,1) = w2fvec6T120(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w2fvec6T120(t,:)'-wplus,1);
    end 
    Return2fvec6T120net(t,1) = (1+Return2fvec6T120gross(t,1))*(1-(cbps/10000)*Turnover)-1;
    %Three-fund
    dc = 0.025;
    c1kvec = 0:dc:1;
    c2kvec = 0:dc:1;
    Uk = zeros(length(c1kvec),length(c2kvec));
    i = 1;
    j = 1;
    for c1 = c1kvec   
        for c2 = c2kvec
            Pk = zeros(T,1);
            for k = 1:kfold
                Xk = Xis;
                Xk(1+(k-1)*T/kfold:k*T/kfold,:) = [];
                Xoosk = Xis(1+(k-1)*T/kfold:k*T/kfold,:);
                muk = mean(Xk)';
                sigmak = cov(Xk,1);
                wgk = inv(sigmak)*onevec./sum(inv(sigmak)*onevec);
                mugk = wgk'*muk;
                wk = (c1/gam)*inv(sigmak)*muk + (c2/gam)*mugk*inv(sigmak)*onevec;
                Pk(1+(k-1)*T/kfold:k*T/kfold,1) = Xoosk*wk;
            end
            Uk(i,j) = mean(Pk) - (gam/2)*var(Pk);
            j = j+1;
        end
        j = 1;
        i = i+1;
    end
    [maxValue,~] = max(Uk(:));
    [index1,index2] = find(Uk == maxValue);
    c1vec6T120(t,1) = c1kvec(1) + (index1-1)*dc;
    c2vec6T120(t,1) = c2kvec(1) + (index2-1)*dc;
    w3fvec6T120(t,:) = (c1vec6T120(t,1)/gam)*invsigmahat*muhat + ...
                      (c2vec6T120(t,1)/gam)*mughat*invsigmahat*onevec;
    Return3fvec6T120gross(t,1) = w3fvec6T120(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
             wplus(i,1) = w3fvec6T120(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w3fvec6T120(t,:)'-wplus,1);
    end 
    Return3fvec6T120net(t,1) = (1+Return3fvec6T120gross(t,1))*(1-(cbps/10000)*Turnover)-1;

end

EU2f1T120gross = 12*(mean(Return2fvec1T120gross)-(gam/2)*var(Return2fvec1T120gross));
EU3f1T120gross = 12*(mean(Return3fvec1T120gross)-(gam/2)*var(Return3fvec1T120gross));
EU2f2T120gross = 12*(mean(Return2fvec2T120gross)-(gam/2)*var(Return2fvec2T120gross));
EU3f2T120gross = 12*(mean(Return3fvec2T120gross)-(gam/2)*var(Return3fvec2T120gross));
EU2f3T120gross = 12*(mean(Return2fvec3T120gross)-(gam/2)*var(Return2fvec3T120gross));
EU3f3T120gross = 12*(mean(Return3fvec3T120gross)-(gam/2)*var(Return3fvec3T120gross));
EU2f4T120gross = 12*(mean(Return2fvec4T120gross)-(gam/2)*var(Return2fvec4T120gross));
EU3f4T120gross = 12*(mean(Return3fvec4T120gross)-(gam/2)*var(Return3fvec4T120gross));
EU2f5T120gross = 12*(mean(Return2fvec5T120gross)-(gam/2)*var(Return2fvec5T120gross));
EU3f5T120gross = 12*(mean(Return3fvec5T120gross)-(gam/2)*var(Return3fvec5T120gross));
EU2f6T120gross = 12*(mean(Return2fvec6T120gross)-(gam/2)*var(Return2fvec6T120gross));
EU3f6T120gross = 12*(mean(Return3fvec6T120gross)-(gam/2)*var(Return3fvec6T120gross));
EU2f1T120net = 12*(mean(Return2fvec1T120net)-(gam/2)*var(Return2fvec1T120net));
EU3f1T120net = 12*(mean(Return3fvec1T120net)-(gam/2)*var(Return3fvec1T120net));
EU2f2T120net = 12*(mean(Return2fvec2T120net)-(gam/2)*var(Return2fvec2T120net));
EU3f2T120net = 12*(mean(Return3fvec2T120net)-(gam/2)*var(Return3fvec2T120net));
EU2f3T120net = 12*(mean(Return2fvec3T120net)-(gam/2)*var(Return2fvec3T120net));
EU3f3T120net = 12*(mean(Return3fvec3T120net)-(gam/2)*var(Return3fvec3T120net));
EU2f4T120net = 12*(mean(Return2fvec4T120net)-(gam/2)*var(Return2fvec4T120net));
EU3f4T120net = 12*(mean(Return3fvec4T120net)-(gam/2)*var(Return3fvec4T120net));
EU2f5T120net = 12*(mean(Return2fvec5T120net)-(gam/2)*var(Return2fvec5T120net));
EU3f5T120net = 12*(mean(Return3fvec5T120net)-(gam/2)*var(Return3fvec5T120net));
EU2f6T120net = 12*(mean(Return2fvec6T120net)-(gam/2)*var(Return2fvec6T120net));
EU3f6T120net = 12*(mean(Return3fvec6T120net)-(gam/2)*var(Return3fvec6T120net));

Table3PanelDTwoFundT120 = ...
[EU2f1T120gross,EU2f2T120gross,EU2f3T120gross,EU2f4T120gross,EU2f5T120gross,...
EU2f6T120gross;EU2f1T120net,EU2f2T120net,EU2f3T120net,EU2f4T120net,EU2f5T120net,...
EU2f6T120net;mean(cvec1T120),mean(cvec2T120),mean(cvec3T120),mean(cvec4T120),...
mean(cvec5T120),mean(cvec6T120)]
Table3PanelDThreeFundT120 = ...
[EU3f1T120gross,EU3f2T120gross,EU3f3T120gross,EU3f4T120gross,EU3f5T120gross,...
EU3f6T120gross;EU3f1T120net,EU3f2T120net,EU3f3T120net,EU3f4T120net,EU3f5T120net,...
EU3f6T120net;mean(c1vec1T120),mean(c1vec2T120),mean(c1vec3T120),mean(c1vec4T120),...
mean(c1vec5T120),mean(c1vec6T120);mean(c2vec1T120),mean(c2vec2T120),mean(c2vec3T120),...
mean(c2vec4T120),mean(c2vec5T120),mean(c2vec6T120)]

%%%%%%%%%%%%%%
%%% T = 240
%%%%%%%%%%%%%%

T = 240;
a1 = (T-N-2)/N;
a2 = a1*(T-N-1)*(T-N-4)/(T-2);
a3 = a2/T;
a4 = ((T-N-1)*(T-N-4)/(T*(T-2)));
rho = N/T;
M = eye(T)-(ones(T,1)*ones(1,T))./T;
NumberWindow = Ttotal-T;

%%% Initialize vectors (c for two-fund, (c1,c2) for three-fund) 
%%% Six calibrations: 1) Normal (exact), 2) t (asymp), 3) t (exact),
%%% 4) Elliptical (asymp), 5) Elliptical (exact), 6) Cross-validation
cvec1T240 = zeros(NumberWindow,1);
c1vec1T240 = zeros(NumberWindow,1);
c2vec1T240 = zeros(NumberWindow,1);
cvec2T240 = zeros(NumberWindow,1);
c1vec2T240 = zeros(NumberWindow,1);
c2vec2T240 = zeros(NumberWindow,1);
cvec3T240 = zeros(NumberWindow,1);
c1vec3T240 = zeros(NumberWindow,1);
c2vec3T240 = zeros(NumberWindow,1);
cvec4T240 = zeros(NumberWindow,1);
c1vec4T240 = zeros(NumberWindow,1);
c2vec4T240 = zeros(NumberWindow,1);
cvec5T240 = zeros(NumberWindow,1);
c1vec5T240 = zeros(NumberWindow,1);
c2vec5T240 = zeros(NumberWindow,1);
cvec6T240 = zeros(NumberWindow,1);
c1vec6T240 = zeros(NumberWindow,1);
c2vec6T240 = zeros(NumberWindow,1);
w2fvec1T240 = zeros(NumberWindow,N);
w3fvec1T240 = zeros(NumberWindow,N);
w2fvec2T240 = zeros(NumberWindow,N);
w3fvec2T240 = zeros(NumberWindow,N);
w2fvec3T240 = zeros(NumberWindow,N);
w3fvec3T240 = zeros(NumberWindow,N);
w2fvec4T240 = zeros(NumberWindow,N);
w3fvec4T240 = zeros(NumberWindow,N);
w2fvec5T240 = zeros(NumberWindow,N);
w3fvec5T240 = zeros(NumberWindow,N);
w2fvec6T240 = zeros(NumberWindow,N);
w3fvec6T240 = zeros(NumberWindow,N);
Return2fvec1T240gross = zeros(NumberWindow,1);
Return3fvec1T240gross = zeros(NumberWindow,1);
Return2fvec2T240gross = zeros(NumberWindow,1);
Return3fvec2T240gross = zeros(NumberWindow,1);
Return2fvec3T240gross = zeros(NumberWindow,1);
Return3fvec3T240gross = zeros(NumberWindow,1);
Return2fvec4T240gross = zeros(NumberWindow,1);
Return3fvec4T240gross = zeros(NumberWindow,1);
Return2fvec5T240gross = zeros(NumberWindow,1);
Return3fvec5T240gross = zeros(NumberWindow,1);
Return2fvec6T240gross = zeros(NumberWindow,1);
Return3fvec6T240gross = zeros(NumberWindow,1);
Return2fvec1T240net = zeros(NumberWindow,1);
Return3fvec1T240net = zeros(NumberWindow,1);
Return2fvec2T240net = zeros(NumberWindow,1);
Return3fvec2T240net = zeros(NumberWindow,1);
Return2fvec3T240net = zeros(NumberWindow,1);
Return3fvec3T240net = zeros(NumberWindow,1);
Return2fvec4T240net = zeros(NumberWindow,1);
Return3fvec4T240net = zeros(NumberWindow,1);
Return2fvec5T240net = zeros(NumberWindow,1);
Return3fvec5T240net = zeros(NumberWindow,1);
Return2fvec6T240net = zeros(NumberWindow,1);
Return3fvec6T240net = zeros(NumberWindow,1);

%%% Rolling window loop
for t = 1:NumberWindow

    if mod(NumberWindow-t,10) == 0
        WindowsLeftT240 = NumberWindow - t
    end
    
    %Inputs
    Xis = R(t:t+T-1,:);
    Xoos = R(t+T,:);
    muhat = mean(Xis)';
    sigmahat = cov(Xis,1);
    invsigmahat = inv(sigmahat);
    theta2hat = muhat'*invsigmahat*muhat;
    theta2hat = ((T-N-2)*theta2hat-N)/T + ((2*(theta2hat)^(N/2)...
                *((1+theta2hat)^(-(T-2)/2)))/(T*beta(N/2,(T-N)/2)*betainc(theta2hat/(1+theta2hat),N/2,(T-N)/2)));
    wghat = invsigmahat*onevec./sum(invsigmahat*onevec);
    mughat = wghat'*muhat;
    varghat = wghat'*sigmahat*wghat;
    psi2hat = muhat'*invsigmahat*muhat - mughat^2/varghat;
    psi2hat = ((T-N-1)*psi2hat-(N-1))/T + ((2*(psi2hat)^((N-1)/2)...
              *((1+psi2hat)^(-(T-2)/2)))/(T*beta((N-1)/2,(T-N+1)/2)*betainc(psi2hat/(1+psi2hat),(N-1)/2,(T-N+1)/2)));
    [~,~,nuhat] = mlstudent(Xis);
    nuhat = min(1000,nuhat);
    tauhat = sum((Xis-muhat').^2,2);
    tauhat = tauhat./mean(tauhat);

    %1) Normal (exact)
    %Two-fund
    cvec1T240(t,1) = a4*theta2hat/(theta2hat+N/T);
    w2fvec1T240(t,:) = (cvec1T240(t,1)/gam)*invsigmahat*muhat;
    Return2fvec1T240gross(t,1) = w2fvec1T240(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
            wplus(i,1) = w2fvec1T240(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w2fvec1T240(t,:)'-wplus,1);
    end 
    Return2fvec1T240net(t,1) = (1+Return2fvec1T240gross(t,1))*(1-(cbps/10000)*Turnover)-1;
    %Three-fund
    c1vec1T240(t,1) = a4*psi2hat/(psi2hat+N/T);
    c2vec1T240(t,1) = a4*(N/T)/(psi2hat+N/T);
    w3fvec1T240(t,:) = (c1vec1T240(t,1)/gam)*invsigmahat*muhat +...
                      (c2vec1T240(t,1)/gam)*mughat*invsigmahat*onevec;
    Return3fvec1T240gross(t,1) = w3fvec1T240(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
            wplus(i,1) = w3fvec1T240(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w3fvec1T240(t,:)'-wplus,1);
    end 
    Return3fvec1T240net(t,1) = (1+Return3fvec1T240gross(t,1))*(1-(cbps/10000)*Turnover)-1;

    %2) t (asymp) 
    y = @(x) (nuhat-2).*rho.*x./(2*(1-rho));
    fun = @(x) y(x).*exp(y(x)).*Expfun(nuhat/2,y(x))-rho;
    eta = fzero(fun,0.5*(1+nuhat/(nuhat-2)));
    varphi = 2*eta^2*(1-rho)/(nuhat-eta*(nuhat-2)); 
    %Two-fund
    cvec2T240(t,1) = (1-rho)^2*theta2hat/((varphi/eta)*theta2hat+rho);
    w2fvec2T240(t,:) = (cvec2T240(t,1)/gam)*invsigmahat*muhat;
    Return2fvec2T240gross(t,1) = w2fvec2T240(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
            wplus(i,1) = w2fvec2T240(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w2fvec2T240(t,:)'-wplus,1);
    end 
    Return2fvec2T240net(t,1) = (1+Return2fvec2T240gross(t,1))*(1-(cbps/10000)*Turnover)-1;
    %Three-fund
    c1vec2T240(t,1) = (1-rho)^2*psi2hat/((varphi/eta)*psi2hat+rho);
    c2vec2T240(t,1) = (1-rho)^2*(eta/varphi)*rho/((varphi/eta)*psi2hat+rho);
    w3fvec2T240(t,:) = (c1vec2T240(t,1)/gam)*invsigmahat*muhat +...
                      (c2vec2T240(t,1)/gam)*mughat*invsigmahat*onevec;
    Return3fvec2T240gross(t,1) = w3fvec2T240(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
            wplus(i,1) = w3fvec2T240(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w3fvec2T240(t,:)'-wplus,1);
    end 
    Return3fvec2T240net(t,1) = (1+Return3fvec2T240gross(t,1))*(1-(cbps/10000)*Turnover)-1;

    %3) t (exact)
    Expk1vec = zeros(nobs,1);
    Expk2vec = zeros(nobs,1);
    Expk3vec = zeros(nobs,1);
    for m = 1:nobs
        lambdahat = diag(sqrt((nuhat-2)./chi2rnd(nuhat,[T,1])));
        B = lambdahat*M*lambdahat;
        oneslambda = lambdahat*ones(T,1);
        Y = randn(T,N);
        Mat1 = inv(Y'*B*Y);
        Mat2 = Mat1^2;
        Expk1vec(m,1) = trace(Mat1);
        Expk2vec(m,1) = trace(Mat2);
        Expk3vec(m,1) = oneslambda'*Y*Mat2*Y'*oneslambda;
    end
    k1 = a1*mean(Expk1vec);
    k2 = a2*mean(Expk2vec);
    k3 = a3*mean(Expk3vec);
    %Two-fund
    cvec3T240(t,1) = a4*k1*theta2hat/(k2*theta2hat+k3*N/T);
    w2fvec3T240(t,:) = (cvec3T240(t,1)/gam)*invsigmahat*muhat;
    Return2fvec3T240gross(t,1) = w2fvec3T240(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
            wplus(i,1) = w2fvec3T240(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w2fvec3T240(t,:)'-wplus,1);
    end 
    Return2fvec3T240net(t,1) = (1+Return2fvec3T240gross(t,1))*(1-(cbps/10000)*Turnover)-1;
    %Three-fund
    c1vec3T240(t,1) = a4*k1*psi2hat/(k2*psi2hat+k3*rho);
    c2vec3T240(t,1) = a4*(k3/k2)*k1*rho/(k2*psi2hat+k3*rho);
    w3fvec3T240(t,:) = (c1vec3T240(t,1) /gam)*invsigmahat*muhat +...
                      (c2vec3T240(t,1)/gam)*mughat*invsigmahat*onevec; 
    Return3fvec3T240gross(t,1) = w3fvec3T240(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
            wplus(i,1) = w3fvec3T240(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w3fvec3T240(t,:)'-wplus,1);
    end 
    Return3fvec3T240net(t,1) = (1+Return3fvec3T240gross(t,1))*(1-(cbps/10000)*Turnover)-1;

    %4) Elliptical (asymp)
    fun = @(x) sum(1./(T-N+x.*tauhat.*N))-1;
    eta = fzero(fun,[1,10]);
    varphi = (1-rho)./(1/eta^2-sum((N.*tauhat.^2)./(T-N+N.*eta.*tauhat).^2));
    %Two-fund
    cvec4T240(t,1) = (1-rho)^2*theta2hat/((varphi/eta)*theta2hat+rho);
    w2fvec4T240(t,:) = (cvec4T240(t,1)/gam)*invsigmahat*muhat;
    Return2fvec4T240gross(t,1) = w2fvec4T240(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
            wplus(i,1) = w2fvec4T240(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w2fvec4T240(t,:)'-wplus,1);
    end 
    Return2fvec4T240net(t,1) = (1+Return2fvec4T240gross(t,1))*(1-(cbps/10000)*Turnover)-1;
    %Three-fund
    c1vec4T240(t,1) = (1-rho)^2*psi2hat/((varphi/eta)*psi2hat+rho);
    c2vec4T240(t,1) = (1-rho)^2*(eta/varphi)*rho/((varphi/eta)*psi2hat+rho);
    w3fvec4T240(t,:) = (c1vec4T240(t,1)/gam)*invsigmahat*muhat +...
                      (c2vec4T240(t,1)/gam)*mughat*invsigmahat*onevec;
    Return3fvec4T240gross(t,1) = w3fvec4T240(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
            wplus(i,1) = w3fvec4T240(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w3fvec4T240(t,:)'-wplus,1);
    end 
    Return3fvec4T240net(t,1) = (1+Return3fvec4T240gross(t,1))*(1-(cbps/10000)*Turnover)-1;
    
    %5) Elliptical (exact)
    lambdahat = diag(sqrt(tauhat));
    B = lambdahat*M*lambdahat;
    oneslambda = lambdahat*ones(T,1);
    Expk1vec = zeros(nobs,1);
    Expk2vec = zeros(nobs,1);
    Expk3vec = zeros(nobs,1);
    for m = 1:nobs
        Y = randn(T,N);
        Mat1 = inv(Y'*B*Y);
        Mat2 = Mat1^2;
        Expk1vec(m,1) = trace(Mat1);
        Expk2vec(m,1) = trace(Mat2);
        Expk3vec(m,1) = oneslambda'*Y*Mat2*Y'*oneslambda;
    end
    k1 = a1*mean(Expk1vec);
    k2 = a2*mean(Expk2vec);
    k3 = a3*mean(Expk3vec);
    %Two-fund
    cvec5T240(t,1) = a4*k1*theta2hat/(k2*theta2hat+k3*N/T);
    w2fvec5T240(t,:) = (cvec5T240(t,1)/gam)*invsigmahat*muhat;
    Return2fvec5T240gross(t,1) = w2fvec5T240(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
            wplus(i,1) = w2fvec5T240(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w2fvec5T240(t,:)'-wplus,1);
    end 
    Return2fvec5T240net(t,1) = (1+Return2fvec5T240gross(t,1))*(1-(cbps/10000)*Turnover)-1;
    %Three-fund
    c1vec5T240(t,1)  = a4*k1*psi2hat/(k2*psi2hat+k3*rho);
    c2vec5T240(t,1) = a4*(k3/k2)*k1*rho/(k2*psi2hat+k3*rho);
    w3fvec5T240(t,:) = (c1vec5T240(t,1) /gam)*invsigmahat*muhat +...
                      (c2vec5T240(t,1)/gam)*mughat*invsigmahat*onevec; 
    Return3fvec5T240gross(t,1) = w3fvec5T240(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
            wplus(i,1) = w3fvec5T240(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w3fvec5T240(t,:)'-wplus,1);
    end 
    Return3fvec5T240net(t,1) = (1+Return3fvec5T240gross(t,1))*(1-(cbps/10000)*Turnover)-1;

    %6) Cross-validation 
    kfold = 5;
    %Two-fund
    dc = 0.01;
    ckvec = 0:dc:1;
    Uk = zeros(length(ckvec),1);
    i = 1;
    for c = ckvec   
        Pk = zeros(T,1);
        for k = 1:kfold
            Xk = Xis;
            Xk(1+(k-1)*T/kfold:k*T/kfold,:) = [];
            Xoosk = Xis(1+(k-1)*T/kfold:k*T/kfold,:);
            muk = mean(Xk)';
            sigmak = cov(Xk,1);
            wk = (c/gam)*inv(sigmak)*muk;
            Pk(1+(k-1)*T/kfold:k*T/kfold,1) = Xoosk*wk;
         end
    Uk(i,1) = mean(Pk) - (gam/2)*var(Pk);
    i = i+1;
    end
    [~,index] = max(Uk);
    cvec6T240(t,1) = (index-1)*dc;
    w2fvec6T240(t,:) = (cvec6T240(t,1)/gam)*invsigmahat*muhat;
    Return2fvec6T240gross(t,1) = w2fvec6T240(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
            wplus(i,1) = w2fvec6T240(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w2fvec6T240(t,:)'-wplus,1);
    end 
    Return2fvec6T240net(t,1) = (1+Return2fvec6T240gross(t,1))*(1-(cbps/10000)*Turnover)-1;
    %Three-fund
    dc = 0.025;
    c1kvec = 0:dc:1;
    c2kvec = 0:dc:1;
    Uk = zeros(length(c1kvec),length(c2kvec));
    i = 1;
    j = 1;
    for c1 = c1kvec   
        for c2 = c2kvec
            Pk = zeros(T,1);
            for k = 1:kfold
                Xk = Xis;
                Xk(1+(k-1)*T/kfold:k*T/kfold,:) = [];
                Xoosk = Xis(1+(k-1)*T/kfold:k*T/kfold,:);
                muk = mean(Xk)';
                sigmak = cov(Xk,1);
                wgk = inv(sigmak)*onevec./sum(inv(sigmak)*onevec);
                mugk = wgk'*muk;
                wk = (c1/gam)*inv(sigmak)*muk + (c2/gam)*mugk*inv(sigmak)*onevec;
                Pk(1+(k-1)*T/kfold:k*T/kfold,1) = Xoosk*wk;
            end
            Uk(i,j) = mean(Pk) - (gam/2)*var(Pk);
            j = j+1;
        end
        j = 1;
        i = i+1;
    end
    [maxValue,~] = max(Uk(:));
    [index1,index2] = find(Uk == maxValue);
    c1vec6T240(t,1) = c1kvec(1) + (index1-1)*dc;
    c2vec6T240(t,1) = c2kvec(1) + (index2-1)*dc;
    w3fvec6T240(t,:) = (c1vec6T240(t,1)/gam)*invsigmahat*muhat + ...
                      (c2vec6T240(t,1)/gam)*mughat*invsigmahat*onevec;
    Return3fvec6T240gross(t,1) = w3fvec6T240(t,:)*Xoos';
    wplus = zeros(N,1);
    if t == 1
       Turnover = 0;  
    else
        for i = 1:N
             wplus(i,1) = w3fvec6T240(t-1,i)*(1+R(T+t-1,i));
        end 
        Turnover = norm(w3fvec6T240(t,:)'-wplus,1);
    end 
    Return3fvec6T240net(t,1) = (1+Return3fvec6T240gross(t,1))*(1-(cbps/10000)*Turnover)-1;

end

EU2f1T240gross = 12*(mean(Return2fvec1T240gross)-(gam/2)*var(Return2fvec1T240gross));
EU3f1T240gross = 12*(mean(Return3fvec1T240gross)-(gam/2)*var(Return3fvec1T240gross));
EU2f2T240gross = 12*(mean(Return2fvec2T240gross)-(gam/2)*var(Return2fvec2T240gross));
EU3f2T240gross = 12*(mean(Return3fvec2T240gross)-(gam/2)*var(Return3fvec2T240gross));
EU2f3T240gross = 12*(mean(Return2fvec3T240gross)-(gam/2)*var(Return2fvec3T240gross));
EU3f3T240gross = 12*(mean(Return3fvec3T240gross)-(gam/2)*var(Return3fvec3T240gross));
EU2f4T240gross = 12*(mean(Return2fvec4T240gross)-(gam/2)*var(Return2fvec4T240gross));
EU3f4T240gross = 12*(mean(Return3fvec4T240gross)-(gam/2)*var(Return3fvec4T240gross));
EU2f5T240gross = 12*(mean(Return2fvec5T240gross)-(gam/2)*var(Return2fvec5T240gross));
EU3f5T240gross = 12*(mean(Return3fvec5T240gross)-(gam/2)*var(Return3fvec5T240gross));
EU2f6T240gross = 12*(mean(Return2fvec6T240gross)-(gam/2)*var(Return2fvec6T240gross));
EU3f6T240gross = 12*(mean(Return3fvec6T240gross)-(gam/2)*var(Return3fvec6T240gross));
EU2f1T240net = 12*(mean(Return2fvec1T240net)-(gam/2)*var(Return2fvec1T240net));
EU3f1T240net = 12*(mean(Return3fvec1T240net)-(gam/2)*var(Return3fvec1T240net));
EU2f2T240net = 12*(mean(Return2fvec2T240net)-(gam/2)*var(Return2fvec2T240net));
EU3f2T240net = 12*(mean(Return3fvec2T240net)-(gam/2)*var(Return3fvec2T240net));
EU2f3T240net = 12*(mean(Return2fvec3T240net)-(gam/2)*var(Return2fvec3T240net));
EU3f3T240net = 12*(mean(Return3fvec3T240net)-(gam/2)*var(Return3fvec3T240net));
EU2f4T240net = 12*(mean(Return2fvec4T240net)-(gam/2)*var(Return2fvec4T240net));
EU3f4T240net = 12*(mean(Return3fvec4T240net)-(gam/2)*var(Return3fvec4T240net));
EU2f5T240net = 12*(mean(Return2fvec5T240net)-(gam/2)*var(Return2fvec5T240net));
EU3f5T240net = 12*(mean(Return3fvec5T240net)-(gam/2)*var(Return3fvec5T240net));
EU2f6T240net = 12*(mean(Return2fvec6T240net)-(gam/2)*var(Return2fvec6T240net));
EU3f6T240net = 12*(mean(Return3fvec6T240net)-(gam/2)*var(Return3fvec6T240net));

Table3PanelDTwoFundT240 = ...
[EU2f1T240gross,EU2f2T240gross,EU2f3T240gross,EU2f4T240gross,EU2f5T240gross,...
EU2f6T240gross;EU2f1T240net,EU2f2T240net,EU2f3T240net,EU2f4T240net,EU2f5T240net,...
EU2f6T240net;mean(cvec1T240),mean(cvec2T240),mean(cvec3T240),mean(cvec4T240),...
mean(cvec5T240),mean(cvec6T240)]
Table3PanelDThreeFundT240 = ...
[EU3f1T240gross,EU3f2T240gross,EU3f3T240gross,EU3f4T240gross,EU3f5T240gross,...
EU3f6T240gross;EU3f1T240net,EU3f2T240net,EU3f3T240net,EU3f4T240net,EU3f5T240net,...
EU3f6T240net;mean(c1vec1T240),mean(c1vec2T240),mean(c1vec3T240),mean(c1vec4T240),...
mean(c1vec5T240),mean(c1vec6T240);mean(c2vec1T240),mean(c2vec2T240),mean(c2vec3T240),...
mean(c2vec4T240),mean(c2vec5T240),mean(c2vec6T240)]

Table3PanelDComplete = [Table3PanelDTwoFundT60;Table3PanelDTwoFundT120;...
    Table3PanelDTwoFundT240;Table3PanelDThreeFundT60;Table3PanelDThreeFundT120;...
    Table3PanelDThreeFundT240]